(library (grid-gen geom)

(export geom:new
        geom:def
        geom:append
)

(import (guile)
        (rnrs)
        (grid-gen main)
        (grid-gen vars)
        (grid-gen math)
        (grid-gen symath)
        (except (srfi srfi-1) map)
        (ice-9 session)
)

(define geom:append (lambda (l v)
  (append l (map (lamcon v) (iota1 (- 3 (length l)))))
))

(define nosimplify-form (lambda (v) (simp-ex (MATH:simplify (norm-form v)))))
;;(define nosimplify-form (lambda (v) (simplify-form v)))

(define geom:vrs (list (lamcon 'X) (lamcon 'Y) (lamcon 'Z)))


(define geom:proj (lambda (gm)
  (lambda (v xi)
    (apply math-cvector
      (map (lambda (x) ;;(MATH:simplify
                         (cons '/ x)
                       ;;)
           )
           (zip (arr-elements ((gm 'proj2) v xi))
                (map (lambda (y) (apply y xi)) (gm 'gks))
           )
      )
    )
  )
))

(define geom:proj2 (lambda (gm)
  (lambda (v xi)
    (apply math-cvector
      (map (lambda (x)
             ;;(map MATH:simplify
               (apply add
                 (map (lambda (y) `(* ,(apply (car y) xi) ,(cadr y)))
                      (zip x (arr-elements v))
                 )
               )
            ;; )
           )
           (gm 'klm)
      )
    )
  )
))

(define geom:bkproj (lambda (gm)
  (lambda (v xi)
    (let* ((apl (lambda (x) (apply x xi)))
           (apm (lambda (v) (map (lambda (x) (apl x)) v)))
           (cf1 (map apl (gm 'gks)))
           (vvs (map (lambda (x) (apply math-cvector (apm x))) (gm 'klm)))
           (cfs (map (lambda (x) (cons '/ x)) (zip (arr-elements v) cf1)))
          )
      ;;(map-ar
      ;;  MATH:simplify
        (apply add
          (map (lambda (x) (apply mul2 x))
               (zip cfs vvs)
          )
        )
      ;;)
    )
  )
))

(define geom:norm (lambda (gm)
  (lambda (d xi)
    (apply math-cvector
           (map (lambda (x)
                  `(/ (apply x xi) (apply (list-ref (gm 'gks) d) xi))
                )
                (list-ref (gm 'klm) d)
           )
    )
  )
))

(define geom:normS (lambda (gm)
  (lambda (d xi)
    (apply math-cvector
           (map (lambda (x) (apply x xi))
                (list-ref (gm 'klm) d)
           )
    )
  )
))

(define geom:new (lambda (f)
  (letrec* 
    ((Q (jakobi-matrix
          (apply math-cvector (apply f geom:vrs))
          (apply math-cvector geom:vrs)
          geom:vrs
        )
     )

     (Q11 Q)
     (K (arr-col Q11 0))
     (L (arr-col Q11 1))
     (M (arr-col Q11 2))
     (KL (sym-vec-mul K L))
     (LM (sym-vec-mul L M))
     (MK (sym-vec-mul M K))
     (Ql `(+ (* ,(list-ref (arr-elements K)  0)  
                ,(list-ref (arr-elements LM) 0)
             )
             (* ,(list-ref (arr-elements K)  1) 
               ,(list-ref (arr-elements LM) 1)
             )
             (* ,(list-ref (arr-elements K)  2) 
                ,(list-ref (arr-elements LM) 2)
             )
          )
     )
     (p-vrs (lambda (v)
       (cond ((is-var? v) (string->symbol (var-name v)))
             ((list? v) (map p-vrs v))
             (#t v)
       )
     ))
     (sfs (lambda (f)
   ;;;    (format #t "\n\nSFS: ~a\n\n" (p-vrs f))
       (let ((f1 (MATH:simplify f)))
         (lambda xl
           (let* ((vl (make-hash-table)))
             (for-each
               (lambda (x) (hash-set! vl (car x) (cadr x)))
               (zip geom:vrs (geom:append xl 0.0))
             )
             (subst-vars f1 vl)
           )
         )
       )
     ))

     (klm (list (map sfs (map nosimplify-form (arr-elements LM)))
                (map sfs (map nosimplify-form (arr-elements MK)))
                (map sfs (map nosimplify-form (arr-elements KL)))
          )
     )
     (KLM (list (map sfs (map nosimplify-form (arr-elements K)))
                (map sfs (map nosimplify-form (arr-elements L)))
                (map sfs (map nosimplify-form (arr-elements M)))
          )
     )
     (Qm (sfs (nosimplify-form Ql)))
     (gks (map sfs (map nosimplify-form (map sym-vec-mod (list LM MK KL)))))
     (basis (map (lambda (vm)
                   (map 
                     (lambda (v)
                       (lambda (xi)
                         (mul2 (apply v xi) `(/ 1.0 ,(apply (cadr vm) xi)))
                       )
                     )
                     (car vm)
                   )
                 )
                 (zip klm gks)
            )
     )

     (self (lambda (cmd)
             (cond ((eq? cmd 'proj)   (geom:proj self))
                   ((eq? cmd 'proj2)  (geom:proj2 self))
                   ((eq? cmd 'bkproj) (geom:bkproj self))
                   ((eq? cmd 'Qm)     Qm)
                   ((eq? cmd 'f)
                    (lambda args
                      (apply f (geom:append args 0.0))
                    )
                   )
                   ((eq? cmd 'gks)    gks)
                   ((eq? cmd 'klm)    klm)
                   ((eq? cmd 'basis)  basis)
                   ((eq? cmd 'norm)   (geom:norm self))
                   ((eq? cmd 'normS)  (geom:normS self))
                   (#t (int-err "Invalig geom cmd: ~a" cmd))
             )
           )
     )
    )
    self
  )
))

(define geom:def (lambda (nd)
  (geom:new (lambda (x y z) `(,x ,y ,z)))
))

)
